Parameter Dependent Problems

63. Parameter Dependent Problems#

We have observed that adding the constraint to the A-block may improve the convergence of saddle-point solvers. If we add the constraint with a large factor $1/\varepsilon$, we may consider the augmented problem as penalty formulation for the constrained minimization problem, and skip the second equation and the Lagrange parameter completely. We call this a parameter dependent problem:

\[ \big(A + \frac{1}{\varepsilon} B^T C^{-1} B \big) u = f \]

Let $V_0 = \operatorname{ker} B$ the null-space of the $B$-matrix.

On the null-space $V_0$ only the $A$-operator is seen. On its complement $V_0^\bot$, the second part $\frac{1}{\varepsilon} B^T C^{-1} B$ dominates for small $\varepsilon$. Unless $V_0$ is trivial, the condition number degenerates with $O(\varepsilon^{-1})$.

The goal is to develop preconditioners for the parameter dependent problems which are robust in $\varepsilon$.

63.1. Dirichlet boundary conditions by penalty:#

Adding the Dirichlet condition $u = u_D$ on $\Gamma_D$:

\[ \int \nabla u \nabla v + \int_{\Gamma_D} \alpha u v = \int f v + \int_{\Gamma_D} \alpha u_D v \]

with large parameter $\alpha = \varepsilon^{-1}$, or depending on the mesh-size $\alpha = h^{-1}$.

The null-space $V_0$ is the finite element sub-space of $H_{0,D}^1 = \{ v \in H^1 : v = 0 \text{ on } \Gamma_D \}$.

63.2. Penalty formulation for the flux:#

Approximating the mixed formulation $\sigma = \nabla u$, and $\operatorname{div} \sigma = -f$ by penalty leads to:

\[ \int \sigma \tau + \frac{1}{\varepsilon} \int \operatorname{div} \sigma \operatorname{div} \tau = \frac{-1}{\varepsilon} \int f \operatorname{div} \tau \]

This is a variational formulation on the Hilbert-space $H(\operatorname{div}) := \{ \tau \in [L_2]^d : \operatorname{div} \tau \in L_2\}$.

The null-space $V_0$ consists of all $\operatorname{div}$-free functions.

63.3. Maxwell equations:#

Very similar is the equation for the magnetic field, a special case of the Maxwell equations. Given a divergence-free current density $j$, search for a divergence-free magnetic field $B$ such that

\[ \operatorname{curl} \mu^{-1} B = j, \]

where $\mu$ is the so called permeability, a kind of conductivity for the magnetic field. Since we search for a $\operatorname{div}$-free $B$, and assuming a simply connected domain, we may introduce a vector-potential $u$ such that $B = \operatorname{curl} u$. Together we have the second order equation

\[ \operatorname{curl} \mu^{-1} \operatorname{curl} u = j. \]

Its weak formulation (skipping discussion of boundary conditions) is

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v = \int j v \qquad \forall \, v \]

There is no unique solution: adding an arbitrary gradient field to a solution is also a solution. One way out is Coloumb - gauging, where we add the constraint $\int u \nabla \psi = 0$, i.e. a weak formulation of $\operatorname{div} u = 0$. Another, very simple method is to regularize with a small $L_2$-term:

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v + \int \varepsilon u v = \int j v \]

Again, we have a parameter dependent problem, with the fictitious regularization parameter $\varepsilon$. The canonical space is

\[ H(\operatorname{curl}) = \{ u \in [L_2]^3 : \operatorname{curl} u \in [L_2]^3 \}, \]

equipped with the norm $$ \| u \|_{H(\operatorname{curl})}^2 = \| u \|_{L_2}^2 + \| \operatorname{curl} u \|_{L_2}^2 $$

The null-space is $$ V_0 = \{ v : \operatorname{curl} v = 0 \} = \nabla H^1 $$

63.4. Penalty formulation for the Stokes equation:#

We approximate the div-free condition by penalty:

Find $u \in [H_{0,D}]^d$ such that

\[ \int \nu \nabla u : \nabla v + \frac{1}{\varepsilon} \int \operatorname{div} u \operatorname{div} v = \int f v \]

It looks similar as the penalty formulation for the flux, but now the A-form is a $H^1$-inner product instead of an $L_2$-inner product. Now, conforming low order methods lead to bad approximations (called locking effect). Robust discretizations are constructed and analyzed by by mixed formulations. If the pressure space consists of discontinuous finite elements, and the $C$ matrix is chosen as $L_2$-inner product, it can be cheaply inverted. We end up with the discretization of the parameter dependent problem. However, the detour over the mixed mixed formulation leads to a projection operator $R_h$:

Find $u_h \in V_h \subset [H_{0,D}]^d$ such that

\[ \int \nu \nabla u_h : \nabla v_h + \frac{1}{\varepsilon} \int R_h \operatorname{div} u_h \operatorname{div} v_h = \int f v_h \]

A stable example is to choose $V_h$ a $P^2$ finite element space. Then $\operatorname{div} u_h$ is a $P^1$ function. Choosing a $P^0$ pressure leads to a projection operator $R_h$ which is the $L_2$-orthogonal projection from $P^1$ onto $P^0$.

64. Robust two-level methods for parameter dependent problems#

Let

\[ A_h^\varepsilon = A + \frac{1}{\varepsilon} B^T C^{-1} B \]

be the discretization matrix of the parameter dependent problem on the finite element space $V_h$. We are going to analyze two-level preconditioners

\[ C_{2L}^{-1} = D_h^{-1} + P \, ( A_H^\varepsilon ) ^{-1} P^T, \]

where $A_H^\varepsilon$ is the discretization matrix on the coarse space $V_H$, $P : V_H \rightarrow V_h$ is the prolongation matrix, and $D_h$ is a local block-Jacobi preconditioner on $V_h$.

64.1. Robust smoothers#

By the ASM - lemma, we can represent the norm generated by the smoother $D_h$ as

\[ \| u \|_{D_h}^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2, \]

where we have the sub-space splitting

\[ V_h = \sum V_i. \]

If $u$ is a function in the finite-element null-space $V_0$, then its norm $\| u \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2$ is independent of the small parameter $\varepsilon$. Only if we can split $u$ into local null-space functions $u_i$, the $D_h$ norm is also independent of $\varepsilon$. Thus, we want a splitting compatible with the null-space:

\[ V_0 = \sum V_i \cap V_0 \]

Example: Maxwell equations

For a compatible discretization with Nedelec finite elements, the discrete null-space consists of discrete gradients:

\[ \{ v_h \in V_h : \operatorname{curl} v_h = 0 \} = \{ \nabla \psi_h : \psi_h \in W_h \subset H^1 \} \]

Let $u_h \in V_0$, and $\psi_h$ such that $\nabla \psi_h = u_h$. Now, decompose $\psi_h = \sum \psi_i$. Then, $\sum \nabla \psi_i$ is a decomposition of $u_h$ into null-space functions $u_i = \nabla \psi_i$. The block-Jacobi must consist of blocks large enough such that each component $\nabla \psi_i$ is contained in one block. If $\psi_i$ is a vertex-hat basis function, then its gradient is contained in the vertex patch. These block-smoothers are called AFW-blocks, due to Arnold, Falk, and Winther: Multigrid in H(div) and H(curl), Numerische Mathematik, 2000.

Lets try:

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))
# mesh.ngmesh.Refine()

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (mesh, clipping=clipping);

fes = HCurl(mesh, order=0)   # one dof per edge
u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx).Assemble()
M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet")).Assemble()

gfu = GridFunction(fes)

pre = Preconditioner(a, "local", block=True)
pre.Update()
from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.46353276292669277     
CG iteration 2, residual = 0.42789858702520756     
CG iteration 3, residual = 0.3232297657475953     
CG iteration 4, residual = 0.254644035863137     
CG iteration 5, residual = 0.18615733922755937     
CG iteration 6, residual = 0.14800181052394998     
CG iteration 7, residual = 0.11949520952076037     
CG iteration 8, residual = 0.09282291331667335     
CG iteration 9, residual = 0.07679038814090586     
CG iteration 10, residual = 0.0638953159561212     
CG iteration 11, residual = 0.05056242146492359     
CG iteration 12, residual = 0.0406053962151821     
CG iteration 13, residual = 0.03430579733378295     
CG iteration 14, residual = 0.030065910434029626     
CG iteration 15, residual = 0.025774846052104745     
CG iteration 16, residual = 0.019839653390976167     
CG iteration 17, residual = 0.013867502737497255     
CG iteration 18, residual = 0.00956747011271997     
CG iteration 19, residual = 0.006601850082304118     
CG iteration 20, residual = 0.004463903340213467     
CG iteration 21, residual = 0.0029382112665803438     
CG iteration 22, residual = 0.002000033787601726     
CG iteration 23, residual = 0.0013982543263260793     
CG iteration 24, residual = 0.0009922984038292016     
CG iteration 25, residual = 0.00078153515076771     
CG iteration 26, residual = 0.0006836420148986603     
CG iteration 27, residual = 0.0006042081173458464     
CG iteration 28, residual = 0.0004933831969725599     
CG iteration 29, residual = 0.00036588144175307274     
CG iteration 30, residual = 0.0002746182002796475     
CG iteration 31, residual = 0.00021606301815491394     
CG iteration 32, residual = 0.00018526566388713884     
CG iteration 33, residual = 0.000167576130198673     
CG iteration 34, residual = 0.00015014932580951798     
CG iteration 35, residual = 0.000126352467517647     
CG iteration 36, residual = 0.00010095216118875678     
CG iteration 37, residual = 7.677150772468239e-05     
CG iteration 38, residual = 5.705147728700163e-05     
CG iteration 39, residual = 4.4009910834998236e-05     
CG iteration 40, residual = 3.641737694358909e-05     
CG iteration 41, residual = 3.212176823485782e-05     
CG iteration 42, residual = 2.7186175188941577e-05     
CG iteration 43, residual = 2.323767518453502e-05     
CG iteration 44, residual = 2.069312434378413e-05     
CG iteration 45, residual = 1.9116273076885043e-05     
CG iteration 46, residual = 1.8718180412921404e-05     
CG iteration 47, residual = 1.8424462724680965e-05     
CG iteration 48, residual = 1.771007517693337e-05     
CG iteration 49, residual = 1.6499366447798923e-05     
CG iteration 50, residual = 1.5042919233742126e-05     
CG iteration 51, residual = 1.4060059508071957e-05     
CG iteration 52, residual = 1.3294831066132331e-05     
CG iteration 53, residual = 1.2250262778532654e-05     
CG iteration 54, residual = 1.0757704058995455e-05     
CG iteration 55, residual = 8.740608113496539e-06     
CG iteration 56, residual = 6.838580733766375e-06     
CG iteration 57, residual = 5.401844236590146e-06     
CG iteration 58, residual = 4.329876676401989e-06     
CG iteration 59, residual = 3.3928176597161652e-06     
CG iteration 60, residual = 2.533031296340977e-06     
CG iteration 61, residual = 1.8791455960249684e-06     
CG iteration 62, residual = 1.512762255204436e-06     
CG iteration 63, residual = 1.4007648574390617e-06     
CG iteration 64, residual = 1.4201528772593683e-06     
CG iteration 65, residual = 1.4845430775582333e-06     
CG iteration 66, residual = 1.5066457507464601e-06     
CG iteration 67, residual = 1.5221995110389638e-06     
CG iteration 68, residual = 1.4801364238907002e-06     
CG iteration 69, residual = 1.41107601815559e-06     
CG iteration 70, residual = 1.253307790798565e-06     
CG iteration 71, residual = 1.0275036694253242e-06     
CG iteration 72, residual = 8.437623855153065e-07     
CG iteration 73, residual = 7.117565782275169e-07     
CG iteration 74, residual = 6.085538137275764e-07     
CG iteration 75, residual = 5.140103061790954e-07     
CG iteration 76, residual = 4.02034778585411e-07     
CG iteration 77, residual = 3.092098916796913e-07     
CG iteration 78, residual = 2.455336702823035e-07     
CG iteration 79, residual = 2.0095472052609272e-07     
CG iteration 80, residual = 1.6554422863396996e-07     
CG iteration 81, residual = 1.3974281461268923e-07     
CG iteration 82, residual = 1.201852101150447e-07     
CG iteration 83, residual = 1.0154715971981071e-07     
CG iteration 84, residual = 8.260511906711917e-08     
CG iteration 85, residual = 6.635817940405876e-08     
CG iteration 86, residual = 5.120263469845712e-08     
CG iteration 87, residual = 4.142946942568193e-08     
CG iteration 88, residual = 3.48569776558864e-08     
CG iteration 89, residual = 2.91715326263099e-08     
CG iteration 90, residual = 2.4620885628215454e-08     
CG iteration 91, residual = 2.114389138424115e-08     
CG iteration 92, residual = 1.8993136309232248e-08     
CG iteration 93, residual = 1.7017039819783714e-08     
CG iteration 94, residual = 1.4333454712749738e-08     
CG iteration 95, residual = 1.1557028423368965e-08     
CG iteration 96, residual = 8.95521894785554e-09     
CG iteration 97, residual = 6.947450640131696e-09     
CG iteration 98, residual = 5.438531508828175e-09     
CG iteration 99, residual = 4.346737344169963e-09     
CG iteration 100, residual = 3.608766183824405e-09     
CG iteration 101, residual = 3.1829807260670588e-09     
CG iteration 102, residual = 2.78436417963747e-09     
CG iteration 103, residual = 2.349635933848687e-09     
CG iteration 104, residual = 2.0666980033295907e-09     
CG iteration 105, residual = 1.8682450158095276e-09     
CG iteration 106, residual = 1.6702589190879603e-09     
CG iteration 107, residual = 1.560887297099862e-09     
CG iteration 108, residual = 1.4141783443008013e-09     
CG iteration 109, residual = 1.180184406983406e-09     
CG iteration 110, residual = 9.680627092761685e-10     
CG iteration 111, residual = 8.362346931082713e-10     
CG iteration 112, residual = 7.562143500985835e-10     
CG iteration 113, residual = 6.519927604250954e-10     
CG iteration 114, residual = 5.603811278094885e-10     
CG iteration 115, residual = 5.069583125986455e-10     
CG iteration 116, residual = 4.926136120245891e-10     
CG iteration 117, residual = 5.022852060108803e-10     
CG iteration 118, residual = 4.873362457890519e-10     
CG iteration 119, residual = 4.1790389403507076e-10     
CG iteration 120, residual = 3.3427726344036365e-10     
CG iteration 121, residual = 2.6056347314262534e-10     
CG iteration 122, residual = 2.0391573266009204e-10     
CG iteration 123, residual = 1.606000240843108e-10     
CG iteration 124, residual = 1.2908545820591397e-10     
CG iteration 125, residual = 1.0736797746566526e-10     
CG iteration 126, residual = 9.326212477559403e-11     
CG iteration 127, residual = 8.43405037487901e-11     
CG iteration 128, residual = 7.529219409679964e-11     
CG iteration 129, residual = 6.44563252530318e-11     
CG iteration 130, residual = 5.515875985155759e-11     
CG iteration 131, residual = 4.642173215709178e-11     
CG iteration 132, residual = 3.969780154434387e-11     
CG iteration 133, residual = 3.510174890456183e-11     
CG iteration 134, residual = 3.133217720518166e-11     
CG iteration 135, residual = 2.75761092369593e-11     
CG iteration 136, residual = 2.354376282472649e-11     
CG iteration 137, residual = 1.8993384753814506e-11     
CG iteration 138, residual = 1.4438615986065556e-11     
CG iteration 139, residual = 1.128483976418334e-11     
CG iteration 140, residual = 9.820233233756264e-12     
CG iteration 141, residual = 9.007363319268733e-12     
CG iteration 142, residual = 8.63350984608252e-12     
CG iteration 143, residual = 8.440755242303782e-12     
CG iteration 144, residual = 7.829576229836768e-12     
CG iteration 145, residual = 6.897243001220731e-12     
CG iteration 146, residual = 5.952563489549341e-12     
CG iteration 147, residual = 5.264312754016778e-12     
CG iteration 148, residual = 4.840281593231537e-12     
CG iteration 149, residual = 4.4923274043033e-12     
CG iteration 150, residual = 4.088488096582208e-12     
CG iteration 151, residual = 3.550560488223952e-12     
CG iteration 152, residual = 2.8329531254937317e-12     
CG iteration 153, residual = 2.1757940510622048e-12     
CG iteration 154, residual = 1.6635119056728308e-12     
CG iteration 155, residual = 1.2634480107340436e-12     
CG iteration 156, residual = 9.361658410851287e-13     
CG iteration 157, residual = 6.973158524308296e-13     
CG iteration 158, residual = 5.230282797657222e-13     
CG iteration 159, residual = 4.113293081154406e-13     

Draw (curl(gfu), mesh, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);

from ngsolve.la import EigenValues_Preconditioner
lam = EigenValues_Preconditioner(a.mat, pre.mat)
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.02952541052575453 5.916054099105063
kappa= 200.37161190170704

We observe that the condition number of the block-Jacobi preconditioner

is robust in $\varepsilon$
depends on the mesh-size

We observe the the condition number of the Jacobi preconditioner

depends on $\varepsilon$
depends on the mesh-size

64.2. Robust coarse-grid correction#

If the penalty term involves a projection, then the coarse-grid null-space $V_{H,0}$ is not a sub-space of the fine-grid null-space $V_{h,0}$. An example is Stokes with penalty:

\[ \| u_H \|_{A_H^\varepsilon}^2 = \| u_H \|_A^2 + \frac{1}{\varepsilon} \| R_H \operatorname{div} u_H \|_{L_2}^2 \]

The operator-norm of the prolongation is $$ \| P \|_{(V_H, A_H^\varepsilon) \rightarrow (V_h, A_h^\varepsilon)} = \sup_{v_H} \frac{ \| P v_H \|_{A_h^\varepsilon}}{\| v_H \|_{A_H^\varepsilon}} $$

To be robust in the parameter, the prolongation operator $P : V_H \rightarrow V_h$ must map the coarse-grid null-space into the fine-grid null-space:

\[ P V_{H,0} \subset V_{h,0} \]

If the null-spaces are nested (as for Maxwell equations), the prolongation is canonical embedding.

For the ASM - analysis we need the existence of stable interpolation operators from fine to coarse: $I_H : V_h \rightarrow V_H$. It must be compatible with null-spaces, i.e. $I_H V_{h,0} \subset V_{H,0}$.

64.3. Two-level analysis for Maxwell equations#

Let $V_h$ be a Nedelec finite element sub-space of $H(\operatorname{curl})$, and consider the bilinear-form

\[ A^\varepsilon(u,v) = \int \operatorname{curl} u \operatorname{curl} v + \varepsilon u v \]

We prove that a two-level method with an AFW block-Jacobi smoother provides an optimal preconditioner.

Let $W_H \subset W_h \subset H^1$ be compatible finite element spaces such that $\nabla W_h = V_{h,0} = \{ v \in V_h : \operatorname{curl} v_h = 0 \}$, and $\nabla W_H = W_{H,0}$.

We need a few technical tools:

Theorem: (Regular Decomposition): There holds

\[ H(\operatorname{curl}) = [H^1]^3 + \nabla H^1 \]

For every function $u$ from $H(\operatorname{curl})$ there exists a decomposition

\[ u = z + \nabla \phi \]

with $z \in [H^1]^3$ and $\phi \in H^1$, and which is stable in the sense

\[\begin{eqnarray*} \| \phi \|_{H^1} & \preceq & \| u \|_{H(\operatorname{curl})} \\ \| z \|_{H^1} & \preceq & \| \operatorname{curl} u \|_{L_2} \\ \end{eqnarray*}\]

On the other hand, every function from $[H^1]^3$, and every gradient from $H^1$ is in $H(\operatorname{curl})$.

Proof: Pasciak, J. E.; Zhao, J.: Overlapping Schwarz methods in H(curl) on polyhedral domains. J. Numer. Math. 10 (2002)

Theorem: (Commuting quasi-interpolation operators) There exist projections $\pi_{V_h} : H(\operatorname{curl}) \rightarrow V_h $ and $\pi_{W_h} : H^1 \rightarrow W_h$ which commute with the gradient:

\[ \pi_{V_h} \nabla = \nabla \pi_{W_h} \]

They are stable in norms and semi-norms, and provide optimal approximation error estimates.

Proof: Schöberl, A multilevel decomposition result in H(curl). in Multigrid, Multilevel and Multiscale Methods, EMG 2005

We are ready to prove that the condition-number of the two-level preconditioner is robust in the mesh-size $h$, as well as the regularization parameter $\varepsilon$. We assume that $H = ch$. We transfer results from the $H^1$-case to the $H(\operatorname{curl})$ case.

In the $H^1$-case, we have shown that the coarse grid interpolation error satisfies

\[ \| u_h - \pi_{W_H} u_h \|_{L_2} \preceq h \| \nabla u_h \| \]

The corresponding lemma for the $H(\operatorname{curl})$ case is:

Lemma: For the coarse-grid interpolation error there exists a decomposition

\[ u_h - \pi_{V_H} u_h = r + \nabla \psi, \]

where $$ \| r \|_{L_2}^2 \preceq h^2 \| \operatorname{curl} u \|^2 $$ and $$ \varepsilon \| \psi \|_{L_2}^2 \preceq h^2 \| u \|_{A^\varepsilon}^2 $$

Proof: Let $u_h \in V_h \subset H(\operatorname{curl})$. There exists a regular decomposition

\[ u_h = z + \nabla \Phi, \]

where

\[ \| z \|_{H^1}^2 \preceq \| \operatorname{curl} u \|_{L_2}^2 \]

and $$ \varepsilon \| \Phi \|_{H^1}^2 \preceq \| u \|_{A^\varepsilon}^2 $$

By the commutativity of the interpolation operator we obtain

\[\begin{eqnarray*} u_h - \pi_{V_H} u_h & = & z - \pi_{V_H} z + \nabla \Phi - \pi_{V_H} \nabla \Phi \\ & = & z - \pi_{V_H} z + \nabla ( \Phi - \pi_{W_H} \Phi ) \end{eqnarray*}\]

By setting

\[ r = z - \pi_{V_H} z \qquad \text{and} \qquad \psi = \Phi - \pi_{W_H} \Phi \]

and using the approximation properties of the interpolation operators we obtain the splitting.

Lemma: For the coarse-grid interpolation error there exists a discrete decomposition

\[ u_h - \pi_{V_H} u_h = r_h + \nabla \psi_h. \]

Proof: Apply the commuting projection operator $\pi_{V_h}$ to the previous lemma to obtain

(64.1)#\[\begin{eqnarray} u_h - \pi_{V_H} u_h & = & \pi_{V_h} ( u_h - \pi_{V_H} u_h ) \\ & = & \pi_{V_h} r + \pi_{V_h} \nabla \psi \\ & = & \pi_{V_h} r + \nabla \pi_{W_h} \psi \end{eqnarray}\]

The discrete decomposition is now $r_h = \pi_{V_h} r$ and $\psi_h = \pi_{W_h} \psi$.

In the $H^1$-case, we can apply the ASM lemma to analyze the Jacobi preconditioner $D$ as follows:

\[ \| u \|_D^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{H^1}^2 \preceq \inf_{u = \sum u_i} h^{-2} \sum \| u_i \|_{L_2}^2 \preceq h^{-2} \| u \|_{L_2}^2 \]

Lemma: For the AFW block-Jacobi preconditioner $D$ there holds

\[ \| u_h \|_D^2 \preceq \inf_{u_h = r_h + \nabla \psi_h} h^{-2} \| r_h \|_{L_2}^2 + \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \]

Proof: Let $u_h = r_h + \nabla \psi_h$ be an arbitrary decomposition. Then the triangle inequality implies

\[ \| u_h \|_D \leq \| r_h \|_D + \| \nabla \psi_h \|_D \]

Now we apply the ASM lemma for both terms:

\[\begin{eqnarray*} \| r_h \|_D^2 & = & \inf_{r = \sum r_i} \sum \| r_i \|_{A^\varepsilon}^2 \\ & = & \inf_{r = \sum r_i} \sum \| \operatorname{curl} r_i \|_{L_2}^2 + \varepsilon \| r_i \|_{L_2}^2 \\ & \preceq & \inf_{r = \sum r_i} h^{-2} \| r_i \|_{L_2}^2 \\ & \approx & h^{-2} \| r \|_{L_2}^2 \end{eqnarray*}\]

For the decomposition of $\nabla \psi$ we use that we can decompose $\psi = \sum \psi_i$ with $\psi_i \in W_i$, where $\nabla W_i \subset V_i$:

\[\begin{eqnarray*} \varepsilon \| \nabla \psi_h \|_D^2 & = & \inf_{\nabla \psi_h = \sum v_i} \sum \| v_i \|_{A^\varepsilon}^2 \\ & \leq & \inf_{\psi_h = \sum \psi_i} \sum \| \nabla \psi_i \|_{A^\varepsilon}^2 \\ & = & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon \| \nabla \psi_i \|_{L_2}^2 \\ & \preceq & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon h^{-2} \| \psi_i \|_{L_2}^2 \\ & \approx & \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \end{eqnarray*}\]

Since $\| u \|_D$ is bounded by any such decomposition of $u_h$, it holds also for the infimum, and the lemma is proven.

Theorem The two-level preconditioner with AFW block-smoothers provide a optimal condition number $\kappa = O(1)$.

Proof: We have to verify that

\[ \inf_{u_h = u_H + \sum u_i} \left\{ \| u_H \|_{A^\varepsilon}^2 + \sum \| u_i \|_{A^\varepsilon}^2 \right\} \preceq \| u_h \|_{A^\varepsilon}^2 \]

We choose $u_H := \pi_{V_H} u_h$. Continuity in energy-norm of the projector bounds the first term. The interpolation rest provides a stable splitting

\[ u_f = u_h - u_H = r_h + \nabla \psi_h, \]

where the components $r_h$ and $\psi_h$ are bounded in $L_2$-norms. This $u_f$ can be further split into local functions as $u_f = \sum u_i$.

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))

fes = HCurl(mesh, order=0, autoupdate=True)
print ("ndof =", fes.ndof)

u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx)
pre = Preconditioner(a, "multigrid", smoother="block")

M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet"))

with TaskManager():
    a.Assemble()
    f.Assemble()
    for l in range(1):  # try also 2
        mesh.Refine()
        print ("ndof =", fes.ndof)
        a.Assemble()
        f.Assemble()

gfu = GridFunction(fes)

ndof = 8800

WARNING: kwarg 'smoother' is an undocumented flags option for class <class 'ngsolve.comp.Preconditioner'>, maybe there is a typo?

ndof = 77372

from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
with TaskManager():
    gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.6907017058419058

CG iteration 2, residual = 0.016142253726008523     
CG iteration 3, residual = 0.000460608017777803     
CG iteration 4, residual = 1.0543115072655466e-05     

CG iteration 5, residual = 3.8277013570456464e-07

CG iteration 6, residual = 3.3711488382831105e-08     
CG iteration 7, residual = 2.4772398927292263e-09     
CG iteration 8, residual = 1.0905496467037034e-10     

CG iteration 9, residual = 1.115369082057783e-11

CG iteration 10, residual = 5.439845265111841e-13

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (curl(gfu), mesh, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);